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Abstract 

Monte Carlo simulation is one of the most important tools in the study of 
i-C ; diffusion processes. For constant diffusion coefficients, an appropriate Gaus- 

^"1 sian distribution of particle's steplengths can generate exact results, when 

Q< compared with integration of the diffusion equation. It is important to no- 

tice that the same method is completely erroneous when applied to non- 
homogeneous diffusion coefficients. A simple alternative, jumping at fixed 
steplegths with appropriate transition probabilities, produces correct results. 
Here, a model for diffusion of calcium ions in the neuromuscular junction of 
the crayfish is used as a test to compare Monte Carlo simulation with fixed 
J^ 4 . and Gaussian steplegth. 

1 Introduction 

Diffusion processes with spatial dependent diffusion coefficient are not uncommon 
in physical, chemical or biological systems. The literature on the subject is vast; 
we mention just a few examples. There is great interest in diffusion processes in 
a tube or channel of varying cross section [IJ[2i[3j|IJ[5j[6j[7J[E]- For a channel 
with cross section depending periodically on the longitudinal direction, it has been 
shown that the dynamics can be described with a one dimensional equation: the 
Fick- Jacobs approximation P, [T], and that this approximation can be improved 
with a spatial dependent diffusion coefficient [2]. In pi)] the authors demonstrated 
that Anderson Localization can be viewed as position-dependent diffusion. In [TTj . 
the antiproton flux from dark matter annihilation-decay is studied using a position- 
dependent diffusion coefficient. In [12] a model of heart muscle is developed using 
a reaction-diffusion system with a spatially varying diffusion coefficient in order to 
investigate arrhythmia. In |13j, experimental results on brownian particles diffusing 
in a confined geometry are reported. The particles are trapped between two nearly 
parallel walls. They measured a spatial dependent diffusion coefficient and a drift in 
the direction of the diffusion coefficient gradient in the absence of an external force 
or concentration gradient. In [TJ] the gating in ion channels is described using a 
position-dependent stochastic diffusion model. 

In [15], a one-dimensional (ID) model for the description of diffusion of Ca 2+ in 
the neuromuscular junction of the crayfish was proposed, in which a fivefold increase 
of the diffusion coefficient takes place over a distance of around 100 nm. 
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Monte Carlo simulation is one of the most important tools in the study of diffu- 
sion processes. It is known that, for a constant diffusion coefficient, the Monte Carlo 
(MC) method with a steplength taken from an appropriate Gaussian distribution 
gives exact results; this means that, on average, the result is the same than the one 
obtained from integration of the diffusion equation. The situation is different for 
non-homogeneous systems. 

The authors of [16] used the model of Ref. [15] to show that the MC method 
with Gaussian steplength and spatial dependent diffusion coefficient (i.e., spatial 
dependent width of the Gaussian distribution) leads to a systematic error, when 
MC averages are compared with the deterministic solution of the diffusion equation. 
To correct this error, they proposed a modification of the Gaussian distribution 
whose main ingredient is to transform it into an asymmetric distribution in order 
to take into account the drift generated by the diffusion coefficient gradient. 

In Ref. [T7], we proposed a different solution. We considered a discrete proba- 
bility distribution for particle jumps at fixed distances to the left and right in a ID 
lattice. The jump length is constant over the whole system, and the jump proba- 
bilities are evaluated in terms of the diffusion coefficient and its gradient through 
simple relations. The method applies to any smooth enough diffusion coefficient. In 
Ref. [IT], we illustrated the method with a constant gradient diffusion coefficient. 
Here, we apply it to the ID model for diffusion of calcium ions in the neuromuscular 
junction of the crayfish [15] and compare the results with the ones obtained from MC 
simulation with Gaussian steplength and from numerical integration of the diffusion 
equation. As already reported in [16], the Gaussian steplength introduces an error 
that even increases with time. The fixed steplength method proposed here gives a 
good agreement when compared with the deterministic results. 

The model is used here as a test for fixed and Gaussian steplength in MC simu- 
lations. Nevertheless, it also has biological interest. Diffusion of calcium ions plays 
an important role in chemical synapses. Chemical synapses are specialized junctions 
through which neurons signal to each other and to non-neuronal cells such as those 
in muscles or glands. The model includes a pulsed input of calcium ions generated 
by an action potential train (the input point is located in the voltage-gated calcium 
channel in the axon terminal). 

The paper is organized as follows. In Sect. [2] we present the general continuous 
and discrete description of the system. In Sect. |3]we describe the particular model 
and the results. In Sect. H]we state our conclusions. 

2 Continuous and discrete descriptions 
2.1 Continuous description 

The connection between the diffusion or Fokker-Planck equation and the stochastic 
motion of Brownian particles is simple when the diffusion coefficient is homogeneous. 
The problems grow when spatial dependence of the diffusion coefficient is taken 
into account. Let us consider a ID system. The simplest Langevin equation for 
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a Brownian particle, in position x, time t, and without a deterministic force (pure 
diffusion) is 

^ = y/2D@jtit), (f) 

where D(x) is the diffusion coefficient and ((t) is the Gaussian white noise assumed 
to be 5-correlated. The corresponding Fokker-Planck equation for the concentration 
of particles c(x, t) is [18] 



dc 
dt 



d_ 
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3D dc 

a ~a~ c + D TT 
ox ox 



(2) 



where a is the interpretation parameter. The simplest interpretations (but not 
the only ones) are the ones of Ito (a = 1), Stratonovich (a = 1/2) and Hanggi- 
Klimontovich (a = 0) [IB] . It is clear from Eq. fl5]) that the interpretation parameter 
is relevant only when D depends on x. 

We consider the standard kinetic equation that is derived from the combination 
of continuity equation and Fick's law, 



dc d 
dt dx 



D 



dc 
dx 



(3) 



that corresponds to a = 0. In this case, the equilibrium solution is homogeneous. 
In a microscopic level this situation can be represented by diffusion in a potential 
where all potential wells are at the same level and the spatial dependence of D is 
originated in variations of barrier's heights (barrier model). This spatial dependence 
produces a drift velocity given by 



dD 

dx 



(4) 



(see Ref. pi]). 



2.2 Discrete description for MC simulation 

A usual choice to simulate diffusion processes is MC method with Gaussian steplength: 
in each MC step time is increased by At and particles jump a distance given by 
y^2D(x)At^, where £ is a Gaussian stochastic variable with mean zero and disper- 
sion 1. As mentioned before, the MC simulation with Gaussian steplength produces 
exact results when the diffusion coefficient is homogeneous. On the other hand, 
also in the homogeneous case, the fixed steplength choice has some important good 
features. It produces exact results in the evaluation of first and second order mo- 
ments [IT] (scaled errors in higher order moments behave as t~ l for t large). And, 
more importantly, it keeps these good features when the diffusion coefficient is non 
homogeneous, in contrast with the Gaussian steplength. 

To implement the fixed steplength simulation, we consider a one dimensional 
discrete lattice with sites separated by potential barriers of different heights. For 
example, Fig. [TJ shows the spatial dependence of the diffusion coefficient of cal- 
cium ions in the neuromuscular junction of the crayfish [15], and, below, we show 
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Figure 1: Top: diffusion coefficient D [nm 2 /ms] of calcium ions against position 
x [/im]. In the fixed steplength simulation the continuous space is replaced by a 
discrete lattice with energy barriers. Bottom: the barriers are shown schematically 
(actual distance between sites is much smaller in the simulations). 
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schematically the corresponding barrier model for diffusion in a lattice. Position of 
site number % is X{ = i Ax, where Ax is the (fixed) lattice spacing. In each MC 
step, a particle in position x, can jump to the right or to the left with probabilities 
p(xi) and q(xi) respectively, and time is increased by a small time step At. The 
transition probabilities per unit time are P(xi) = p(xi)/At and Q(xi) = q{xj)/At. 
Using the corresponding Master equation for these transition processes, we can ar- 
rive to a discretized version of the diffusion equation (jHJ), from which we can write 
the relations between transition probabilities and diffusion coefficient (see Eqs. (7) 
and (8) of Ref. [HI): 

„. , D(xi) 1 3D, . 

The average drift velocity can be obtained from these equations: 

dD 

v = [P(xi) - Q(xi)} Ax = -^-(xi), (7) 

ox 

that is in agreement with Eq. (j4|). The asymmetry between right and left jumps 
produces the drift current. The main problem with the Gaussian distribution of 
steplengths is its symmetry, that precludes the appearance of the drift current. 



3 Model and results 

The model for the diffusion of clacium ions in the neuromuscular junction of the 
cryfish introduced in [15], and used in [16], was proposed in order to represent some 
experimental results for which it was necessary to assume a fivefold increase in the 
diffusion coefficient over a distance of the order of 100 nm. The diffusion coefficient 
is 

D(x)=D[l-0.Bu(x)] (8) 

u(x) = - {tanh [A (6 - x)] + 1} , (9) 

where D = 4 x 10~ 3 /im 2 /ms, A = 35 /im -1 and b = 0.2 /im (see Fig. [T]). The 
point where ions are injected in the system is x = 0, this is where the voltage-gated 
calcium channel is located. An action potential train of frequency 100 Hz opens 
the channel and, for each pulse, 8000 particles of calcium ions enter the system. 
Particles are released during a span of 1.2 ms. The boundary condition at x = is 
reflecting. 

For the MC simulation with Gaussian or fixed steplength we used At = 0.00015 
ms. During the 1.2 ms that lasts a pulse of the action potential, one particle enters 
the system in each MC step, so that 8000 particles enter for each pulse. Results 
were averaged over 50 samples. 

For fixed steplength we used a lattice where position Xj is given by iAx with 
% = 0,...,10 4 and Ax = 0.002 /im (the lattice is large enough to avoid boundary 



5 



3xl0 5 1 r 




t t 

Figure 2: Concentration of particles c (number of particles per /im) against time t 
[ms] evaluated at different distances x to the right of the input site. From left to right 
and from top to bottom: x=0.02, 0.06, 0.3 and 0.55 jum. Squares: fixed steplength; 
continuous curve (over squares): numerical integration of diffusion equation; plus 
signs: Gaussian steplength. 

effects). The value of Ax should be much smaller than any characteristic length of 
the system (in our case, around 0.1 /im) and, combined with At, they should satisfy 
p(xi) + q(xi) < 1 Vz. For cases where jumping probabilities are very small in parts 
of the system, it is convenient to combine the fixed steplength with the exponential 
distribution of time steps given by kinetic MC to reduce computational time [T9] . 

We assume that the system is exactly described by the diffusion equation ([3]). 
The accuracy degree of the MC simulation results is determined comparing them 
with numerical integration of ([3]). Numerical integration was performed using the 
FTCS method [20J. The initial condition is c(x,0) = 0, and the discretization for 
numerical integration is also Ax = 0.002 fim. 

Fig. H] shows the particle concentration as function of time evaluated at different 
distances x to the right of the input site (x = 0.02, 0.06, 0.3 and 0.55 /xm). Numerical 
results with fixed steplength coincide with the results of integration of the differential 
equation for all times and distances. Gaussian steplength results do not coincide, 
and the error increases with time. The error is positive for short distances and 
negative for large distances. This fact can be understood analyzing the effect of the 
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Figure 3: Concentration c (number of particles per /xm) against position x \pm] for 
different times t [ms]. From left to right and from top to bottom: t — 5, 15, 45 
and 55 ms. Squares: fixed steplength; continuous curve (over squares): numerical 
integration of diffusion equation; plus signs: Gaussian steplength. 
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drift current. The spatial dependence of D (Fig. CD) induces a positive drift velocity, 
see Eq. (BJ. For Gaussian steplength the drift is absent, and this produces an excess 
of particles close to the input site (x small) and a lack of particles far from it to the 
right (x large). 

The difference between fixed and Gaussian steplength, and the absence of drift 
in the last one, is clearer in Fig. [3J where particle concentration against position for 
different times t (t = 5, 15, 45 and 55 ms) is plotted. The times were chosen in 
the middle of consecutive pulses. As before, fixed steplength results coincide with 
the deterministic result of the diffusion equation, while Gaussian steplength results 
do not coincide. The curve of fixed steplength results tends to favour larger values 
of c for large x because of the positive drift in the region around x = 0.2 /im. For 
example, in the plot for t = 45 ms or 55 ms, we can clearly see that the Gaussian 
steplength curve has an excess (lack) of particles for x small (large), as mentioned 
before. 

The effect of the absence of drift in the Gaussian steplength results is also present 
in the evaluation of the average position (x) against time, as shown in Fig. HI It can 
be seen that the results of (x) for fixed steplength and numerical integration coincide. 
But, for Gaussian steplength, (x) is smaller for all t. In both cases, (x) increses with 
time due to the continuous input of particles and the reflecting boundary condition 
at x = 0. The difference between Gaussian and fixed steplength is due to the drift 
velocity. Particle input is actually pulsed, but pulses are not present in this plot 
because time was chosen in the middle of consecutive pulses. 

4 Conclusions 

We have used a model for diffusion of calcium ions in the neuromuscular junction of 
the crayfish to compare MC simulation results with fixed and Gaussian steplength. 
The system is assumed to be exactly described by the diffusion equation. Com- 
parison of the MC simulation results with numerical integration of the diffusion 
equation shows a time increasing error for the Gaussian steplength choice, while 
fixed steplength coincide with numerical integration. Results of the particle den- 
sity at different times and at different positions for Gaussian steplength show errors 
that can be understood in terms of the symmetric shape of the Gaussian distri- 
bution, that can not produce the drift generated by the diffusion gradient. This 
feature makes the Gaussian steplength completely inadequate for the simulation of 
non- homogeneous systems. 

In a previous work [17J, we have shown that the fixed steplength method exactly 
reproduces first and second order moments in homogeneous systems, and that the 
scaled error in higher order moments behaves as t _1 for t large. 

The method is proposed as a simpler alternative to the methods of Refs. [TH [21] 
to correct the errors introduced by the Gaussian steplength in non homogeneous 
systems. 
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Figure 4: Average position (x) [//m] against time [ms]. Squares: fixed steplength; 
stars (over squares): numerical integration of diffusion equation; plus signs: Gaus- 
sian steplength. 
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